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Complex network theory aims to model and analyze complex systems that consist of mul- 
tiple and interdependent components. Among all studies on complex networks, topological 
structure analysis is of the most fundamental importance, as it represents a natural route to 
understand the dynamics, as well as to synthesize or optimize the functions, of networks. A 
broad spectrum of network structural patterns have been respectively reported in the past 
decade, such as communities, multipartites, hubs, authorities, outliers, bow ties, and others. 
Here, we show that most individual real- world networks demonstrate multiplex structures. 
That is, a multitude of known or even unknown (hidden) patterns can simultaneously situate 
in the same network, and moreover they may be overlapped and nested with each other to 
collaboratively form a heterogeneous, nested or hierarchical organization, in which different 
connective phenomena can be observed at different granular levels. In addition, we show 
that the multiplex structures hidden in exploratory networks can be well defined as well 
as effectively recognized within an unified framework consisting of a set of proposed con- 
cepts, models, and algorithms. Our findings provide a strong evidence that most real-world 
complex systems are driven by a combination of heterogeneous mechanisms that may col- 
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laboratively shape their ubiquitous multiplex structures as we observe currently. This work 
also contributes a mathematical tool for analyzing different sources of networks from a new 
perspective of unveiling multiplex structures, which will be beneficial to multiple disciplines 
including sociology, economics and computer science. 

1 Introduction 

Complex network analysis provides a novel approach to examining how networked systems in 
nature are originated and evolving according to what basic principles, and moreover armed with 
such discovered principles, constructing efficient, robust as well as flexible man-made networked 
systems under different constraints. Among all studies about complex networks, structure analysis 
is the most fundamental, and the ability to discover and visualize the underlying structure of a 
real- world network in question will be greatly helpful for both topological and dynamic analysis 
applied to it 4 So far, scientists have uncovered a multitude of structural patterns ubiquitously 
existing in social, biological, ecological or technological networks; they may be microscopic such 
as motifs 2 , mesoscopic such as modularitieP or macroscopic such as small world ^ and scale- 
free phenomena 5 . These structural patterns observed at different granular levels may collectively 
unveil the secrets hidden in complex networked systems. All these works have greatly triggered 
the common interesting as well as boosted the progress of exploring complex networks in both 
scientific and engineering domains. However, the topological structure analysis of complex net- 
works, even restricted to the mesoscopic level, remains one of the major challenges in network 
theory mainly because most of the real- world networks are usually resulted from a combination 
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of heterogeneous mechanisms which may collaboratively shape their non-trivial structures. More 
specifically one can raise the following issues. 

Above all, beyond modularity the most extensively studied at the mesoscopic level a wide 
spectrum of structural patterns have been reported in the literature including bipartites or more 
generally multipartites hubs , authorities and outliers bow ties or others. Moreover, 
these miscellaneous patterns may simultaneously coexist in the same networks, or they may even 
overlap with each other such as the overlaps between communities 16 . Fig. [I] shows an illustrative 
example, in which a social network encoding the co-appearances of 77 characters in the novel "Les 
Miserables" is visualized in terms of both network and matrix representations. One would observe 
two hubs and a number of outliers coexisting with four well-formed communities. The two hubs, 
corresponding to Valjean and Javert, are the main characters of the novel, and their links are across 
all other clusters by connecting about 48% of the overall characters. It indicates that the two roles 
interacting with different characters in different chapters are the two main clues going through the 
whole story. Four detected communities can be seen as four relatively independent social cliques. 
As an example, we go into details of group 4 that is almost separated from the rest of the network. 
Interestingly, this small social clique consists of 4 parisian students Tholomyes, Listolier, Fameuil 
and Blacheville with their respective lovers Fantine, Dahlia, Zephine and Favourite. Group 5 
consists of 39 outliers, which connect to either two hubs or one of four communities by only a few 
links. Correspondingly, these outliers are the supporting roles of this novel. Besides this example, 
more complex structural patterns in real- world networks will be demonstrated and discussed in the 
rest of this article. 
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Secondly, such multiple structural patterns may be nested. That is, real- world networks can 
contain hierarchical organizations with heterogeneous patterns at different levels. In the literature, 
hierarchical structures are usually studied in a homogeneous way, in which patterns observed in 
each layer of hierarchies show homophily in terms of either fractal property ^ or modularity in 
more general cases G2H3 A recent study reveals that the patterns demonstrated in each layer of a 
dendrogram can be assortative (or modular) or disassortative (or bipartite)^. The ability to observe 
heterogeneous patterns at different levels brings new clues for understanding the dynamics of real- 
world complex networks. 

Due to the above two reasons, for an exploratory network about which one often knows 
little, it is very hard to know what specific structural patterns can be expected and then be obtained 
by what specific tools. Biased results will be caused if an inappropriate tool is chosen; even if 
we know something about it beforehand, it is still difficult for a tool exclusively designed for 
exploring very specific patterns, say modularity, to satisfactorily uncover a multiple of coexisting 
patterns possibly overlapped or even nested with each other (we call them multiplex structures in 
this work) from networks. To the best our knowledge, there have been no studies in the literature 
being able to address both of the above issues. On the other hand, human beings have the capability 
of simultaneously discovering multiplex and significant structural patterns for various objectives. 
It has been believed that this kind of capability is an important form of human cognitive and 
intelligent functions^. 

Back to the matrix as shown in the Fig. [TJ one would observe an intuitive phenomenon: Any 
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significant pattern contained in the underlying structures of the network can be statistically high- 
lighted by a group of homogenous individuals with identical or quite similar connectivity profiles. 
For instance, individuals from the same communities prefer to intensively interact with each other 
but rarely interact with the rest; hubs would prefer to connect many individuals from different parts 
of the whole network, whereas all outliers tend to seldom play with others by emitting only a few 
connections. Based on this naive observation, if one can group the majority of individuals into 
reasonable clusters according to their connectivity profiles, the coexisting structural patterns can 
be unveiled by further inferring the linkage among clusters. In this way, the first issue listed above 
can be promisingly solved. 

The idea of grouping nodes into equivalent clusters in terms of their connection patterns is 
quite similar with the philosophy of the blockmodeling first proposed by Lorrain and White 23 in 
which nodes with structural equivalence (defined in terms of local neighborhood configurations) 
or more generally regular equivalence ^ (defined in terms of globally physical connections to all 
others) or more softly stochastic equivalence^^ (defined in terms of linking probabilities between 
groups) will be grouped into the same blocks. 

Based on the same idea, a very related study has been proposed recently by Newman and 
Leicht, which first (to our knowledge) shows the motivation to detect unpredefined structural pat- 
terns from exploratory networks 27 . From the perspective of machine leaning, their method can be 
seen as a version of naive Bayes algorithm applied to networks, whose objective is to group nodes 
with similar connection features into a predefined number of clusters. Although their work only 



5 



shows the ability to determine whether an exploratory network is assortative or disassortative by 
manually analyzing the obtained clusters, it has provided one good proof supporting that the idea 
of grouping nodes into equivalent clusters can be an initial step of the whole process aiming to 
unveil multiplex structures from networks. 

In this work, we will propose a novel model by introducing the concept of granularity into 
stochastic connection profiles in order to properly model multiplex structures, and then show how 
the task of recognizing multiplex structures can be reduced to a simplified version of the iso- 
morphism subgraph matching problem. To test our ideas and strategies proposed here, different 
sources of networks have been analyzed. It is encouraging that our methods show a good perfor- 
mance, capable of uncovering multiplex structures from the tested networks in a fully automatic 
way. 



2 The Model 



Granular couplings We will model the connection profiles of nodes in terms of probabilities 
instead of physical connectivity. In this way, it is expected not only to find out multiplex structures 
but also to provide an explicit probabilistic interpretation for these findings within the Bayesian 
framework. We term such probabilities as couplings in that they are not just the mathematical 
measures subjectively defined for modeling or computing, but they do exist in many real- world 
systems, encoding different physical meanings such as social preferences in societies, predation 
habits in ecosystems, co-expression regularities in gene networks or co-occurrence likelihood of 
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words in languages, which will be valuable to predict their situated systems. 

Here the notation of granularity should be interpreted in terms of the resolution and preci- 
sion. On one hand, in our model we define two kinds of couplings with different resolutions, i.e., 
node couplings and block couplings. Formally, we define node feedforward-coupling matrix P nXn 
and node feedback-coupling matrix Q nX m where pij and respectively denote the probabilities 
that node i expects to couple with or to be coupled by node j. In the cases of indirected networks 
we have P = Q (see SI for proofs). We assume nodes will independently couple with others 
regulated by such couplings. Nodes with similar feedforward- as well as feedback- coupling dis- 
tributions will be clustered into the same blocks. In terms of matrices, homogeneous feedforward- 
and feedback-couplings guarantee homogeneous row and column connection profiles, respectively. 
Correspondingly, given the block number K, we define block feedforward-coupling matrix §kxk 
and block feedback-coupling matrix ^kxk, where (j) pq and ip pq respectively denote the probabilities 
that block C p expects to couple with or to be coupled by block C q . We will show later that block 
couplings can be inferred from node couplings and vice versa. Node couplings with a fine granu- 
larity are used to model networks in order to capture their local information as much as possible; 
while block couplings with a coarser granularity are used to define and recognize global structural 
patterns by intentionally neglecting trivial details. On the other hand, in the nested patterns, node 
couplings and block couplings in different hierarchies will have different precisions in order to 
properly abstract and construct hierarchical organizations. In our model, the couplings on higher 
layers are the approximations of the related ones on the lower layers. Therefore, as the layers 
moving from the bottom (corresponding to the original networks) to the top (corresponding to the 
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finally reduced networks) of the nested organizations, the precision of node or block couplings will 
gradually degenerate. 

Defining multiplex structures The main steps of our strategies for detecting multiplex struc- 
tures from networks can be stated as follows: 1) simultaneously estimating all kinds of couplings 
mentioned above and clustering all nodes into nested blocks with a proposed granular blocking 
algorithm; and 2) in each layer of the nested blocks, recognizing structural patterns by matching 
predefined isomorphism subgraphs from a reduced blocking model in which trivial couplings are 
neglected, as illustrated in Fig. [TJi. Multiplex structures can be defined in terms of blocks and 
their couplings. Fig. [2] shows a schematic illustration by means of some conceptual networks. By 
referring to them, we give following definitions. 

A community is a self-coupled block. An authority is a self-coupled block which is coupled 
by a number of other blocks. A hub is a self-coupled block which couples with a number of other 
blocks. An outlier is a block without self-coupling which is coupled by a hub or couples with an 
authority. A bow-tie is a subgraph consisting of a block b and two sets of blocks Si(b) and S r (b), 
which satisfy with: 1) b is coupled by and couples with the blocks of Si(p) and S r (p), respectively; 
2) the intersection of Si(p) and S r (p) is empty or {6}; and 3) there are no couplings between Si(p) 
and S r (p); A multipartite is a subgraph consisting of a set of blocks without self-couplings which 
reciprocally couple with each other. As a special case of multipartite, a bipartite is a subgraph 
consisting of two blocks without self-loop couplings which unilaterally or bilaterally coupled with 
each other. A hierarchical organization is a set of nested blocks, in which block subgraphs in 



8 



lower layers are directly or indirectly nested in the bigger blocks on higher layers. 

The above definitions imply that there may exist overlaps between different patterns in the 
sense that the same blocks can be simultaneously involved in a multitude of subgraphs. For exam- 
ple, a block which is determined as a community can be also a hub, a authority or the core of one 
bow tie. Moreover, beyond the predefined patterns, users are allowed to define novel even more 
complex patterns by designing new subgraphs of blocks, which can be identified by matching their 
isomorphism counterparts from blocking models. 

Granular blocking model Let TV = (V, E) be a directed and binary network, where V(N) de- 
notes the set of nodes and E(N) denotes the set of directed links. In the case of undirected network, 
we suppose there are two direct links between each pair of nodes. Let A nxn be the adjacency matrix 
of N, where n is the number of nodes. 

Suppose all nodes of TV are divided into L(l < L < n) blocks, denotes by B nxL , where 
bu = 1 if node i is in the block /, otherwise bu = 0. When each block is considered to be 
inseparable, the granularity of network TV can be measured by the average size of blocks g = n/L. 
As g increasing from 1 to n, the granularity of TV degenerates from the finest to the coarsest. Let 
B g denote the blocking model with a granularity g, and we expect to cluster all its blocks into a 
reasonable number of clusters so that the nodes of blocks within the same cluster will demonstrate 
homogeneous coupling distributions. Let matrix Z LxK (l < K < L) denote such clusters, where 
K is the cluster number and z ik = 1 if block I is labeled by cluster h, otherwise z\ k = 0. Since the 
coupling distributions of nodes within the same clusters are expected to be homogeneous, one can 

9 



characterize such distributions for each cluster instead of for each node. Given Z, define 9^ xn , 
where 6 k j denotes the probability that any node out of cluster k expects to couple with node Vj ; 
define A Kxn , where 5 k j denotes the probability that any node out of cluster k expects to be coupled 
by node Vj\ define Q = • • • , ujk) t , where u k denotes the prior probability that a randomly 
selected node will belong to the cluster C k . It is easy to show P = B g ZQ and Q = B g ZA (see 
SI). 

Let X = (K, Z, 0, A, Q) be a model with respect to TV and B g . We expect to select an 
optimal X from its hypothesis space to properly fit as well as to precisely predict the behaviors 
of TV under B g in terms of node couplings characterized by it. According to the MAP principle 
(maximum a posteriori), the optimal X for a given network TV under B g will be the one with 
the maximum posterior probability. Moreover, we have: P(X\N J B g ) oc P(N\X J B g )P{X\B g ), 
where P(X\N, B g ), P(N\X, B g ) and P{X\B g ) denote the posteriori of X given N and B g , the 
likelihood of TV given X and B g and the priori of X given B g , respectively. 

3 Learning Methods 

Likelihood maximization We first consider the simplest case by assuming all the prioris of X 
given B g are equal. In this case, we have: P(X\N, B g ) oc P(N\X,B g ). Let L(N\X, B g ) = 
\nP(N\X, B g ), we have (see SI): 
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L K 

L(N\X,B g ) = J2Y1 ln E n M(^%)/fe a i>fc (!) 

Z=l 6^/0 fc=l 

where /(#, y) = — 



Let L(iV, 5 P ) be the log-likelihood of the joint distribution of TV and Z given X and 
£?3, we have (see SI): 



L K n 

L(N, Z\X, B g ) = J2J2Y, M^QnfVkj, (Hi) + \nf(5 kj , a 3l )) + In a;*) (2) 

1=1 b u ^0 k=l j=l 



Considering the expectation of L(N, Z\X, B g ) on Z, we have: 



L K n 

E[L(N,Z\X,B g )] = Y,Y, X^i*E0n^ (3) 

1=1 ba^O k=l j=l 



where E[zi k ] = 7^ = P(y = k\b = /, X, B g ), i.e., the probability that block I will be labeled as 
cluster k given X and B g . 



Let J = E[L(7V, Z\X, B g )\ + Af^^c^ — 1) be a Lagrangian function constructed for 
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maximizing E[L(N, Z\X, B g )] with a constraint J2k=i u k = 1> we have: 





dO k j 



dS kj 



dJ 

duj k 

dJ_ 

d\ 



= 



U kj — spL v — — 

c __ J2iLiJ2b u ^o a j^ik 

°kj — v — ~ 

2^ = 1 2^b u ^0 'Ik 

^fc " 



(4) 



According to the Bayesian theorem,we have (see SI): 



(5) 



Using the similar treatment as proposed by Dempster and Laird ^ we can prove that a local 
optimum of Eq(l0]will be guaranteed by recursively calculating Eqsj4]and[T2](see SI). The time 
complexity of this iterative computing process is 0(In 2 K), where I is the iterations required for 
convergence, which is usually quite small. An approximate but much faster version with a time 
O(ILnK) is given in the SI. 



Priori approximation Without considering the priori of the model, the above-proposed likelihood 
maximization algorithm will be cursed by the overfitting problem. That is, L(N\X, B g ) will mono- 
tonically increase as K approaching to L. In this section, we will discuss how to approximately 
estimate the prior P(X\B g ) by means of the information theory. 

Note that 1 < K < L = n/g, which implies that the coarser granularity the smaller K. It 
will be shown in the following that a smaller K will indicate a less complexity of X. So, we have: 
a coarser granularity prefer simpler models, which can be mathematically written as P(X\B g ) = 
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r](X) 9 , where the function rj(X) measures the complexity of X in terms of its parameters. In this 
work, we set rj(X) = P{X\B\) = P(X) 9 where P(X) is the priori of X in the hypothesis space 
in which K can be freely valued from 1 to n. According to Shannon and Weaver 29 , ln(l/P(X)) 
is the minimum description length of X with a prior P(X) in its model space. Let OC{X) denote 
the optimal coding schema for X, and let L Q c{X) be the minimum description length of X under 
this schema. We have: — \nP(X\B g ) = — g\nP(X) = gL c{X). Now, to estimate the prior 
P{X\B g ) is to design a good coding (or compressing) schema as close to OC as possible. 

In terms of the parameter of X, i.e.,©, A, Z and Q, we have (see SI): 



$ = GB g ZD~\ $ = ABgZD' 1 (6) 

where D = diag(nQ). 

Note that matrices Z and B g can be compressed into a map y, where y(i) = k if the entry 
(z, k) of B g Z is equal to one. Given y, $ and ^, node couplings and g i<7 - can be measured by: 



Pij = </>y(i),y(j)i % = ^y(i),y(j) ( 7 ) 

EqjTjsays that, all node couplings can be approximately characterized by y, $ and ^. In 
other word, the compressing schema close to OC(X) we have found out is: 

OC(X) = (AT, $XXK, *KXX, J/nx2) (8) 
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Now, we have, 



LggpQ = 1 x (- In ±) + 2K\- In + 2n(- In 
= 2K 2 lnK 2 + 2nln2n 

Moreover, we have: 

arg max x P(X|A^, 5 p ) 

= argmin x (-lnP(A^|X,^) -lnP(X|^)) 

= arg min x (- lnP(N\X, B g ) + gL^{X)) (9) 
= arg min x {-L(N\X, B g ) + #(2if 2 In if 2 + 2nln2n)) 
= arg min x (-L(AT|X,^) + 2gK 2 InK 2 ) 

Eqj9] tries to seek a good tradeoff between the accuracy of model (or the precision of fitting 
observed data) measured by the likelihood of a network, and the complexity of a model (or the 
generalization ability to predict new data) measured by its optimal coding length. 

Model selection For a given K, the penalty term is a constant, and thus to maximize P(X\N, B g ) 
is to maximize L(N\X, B g ). In our algorithm, K will be systematically checked from 1 to L, 
and the model with the minimum sum of negative likelihood and penalty will be returned as the 
optimal one. In the landscape of K and P(X\N,B g ) 9 a well-like curve will be shaped during the 
whole search process (see SI). So, in practice, rather than mechanically checking K for exact L 
times, an ongoing searching can be stopped after it has safely passed the well bottom. By means 
of this greedy strategy, the efficiency of our algorithm would be greatly improved. The complete 
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algorithm is given in the SI. 

Hierarchy construction Assume we have constructed h layers, in which the i — th layer corre- 
sponds to a blocking model characterized by B 9i . Now, we want to construct the (h + 1) — th 
layer by selecting an X with a maximum posterior given a set of blocking models on differ- 
ent layers. We have shown that (see SI): P(X\N, B gi , • • • , B 9h ) oc P(N\X, B 9h )P(X\B 9h ) oc 
P(N\X, B gh )P(X) 9h . This Markov-like process indicates that the new model to be selected for 
layer h is only based on the information of the layer h — 1. So, for a given network, its hierarchical 
organization can be incrementally constructed as follows. 

Firstly, we construct the first layer by taking each node as one block, and cluster it into K\ 
clusters by selecting an model with a maximum P{X\N,B\). Next, we form B n / Kl by capsuling 
each cluster on the first layer as one block, and cluster these blocks into K 2 clusters by selecting an 
new model with a maximum P(X\N J B n / Kl ) 9 which forms the second layer. We repeat the same 
procedure until it converges (the cluster number obtained keeps fixed). In this way, the number of 
layers of a hierarchical organization will be automatically determined. The above procedure can 
be seen as a semi- supervised learning process; in the cases that the granularity is more than one, 
we have already known a priori that which nodes will be definitely within the same clusters. As the 
layer in the constructed hierarchy increases, the homogeneity of the nodes within the same cluster 
in terms of their connection profiles keeps degenerating, and correspondingly more global patterns 
are allowed to be observed by tolerating such increasing diversities. 
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Isomorphism subgraph matching Based on the obtained blocks in the level In + 1 of the hierar- 
chy, all potential patterns hidden in the level h can be revealed through an isomorphism subgraph 
matching procedure, whose input is the block feedforward-coupling matrix $. First we construct a 
reduced blocking model by taking each block as one node, and for each pair of blocks p and q we 
insert a link from p to q if cj) pq is above a threshold computed based on $ (see SI). Then for each 
block, we pick up the matched isomorphism subgraphs it will be involved in, and put them into 
categorized reservoirs labeled by different patterns. During this procedure, the subgraph to be put 
into a reservoir will be discarded if it is a subgraph of an existing one, as illustrated by Fig|3]± 

4 Applications 

Exploring the cash flow patterns of the world trade system The discovered multiplex structures 
as well as their granular couplings can be used to understand some dynamics of networks. Here 
we give one example to show how cash possibly flow through a world trade net. Fig. [3^ shows 
a directed network encoding the trade relationship among eighty countries in the world in 1994, 
which was originally constructed by Nooy^ based on the commodity trade statistics published by 
the United Nations. Nodes denote countries, and each arc i —> j denotes the country i imported 
high technology products or heavy manufactures from the country j. Analogous to the "structural 
world system positions" initially suggested by Smith and White based on their analysis of the 
international trade from 1965 to 1980^, the eight countries in 1994 were categorized into three 
classes according to their economic roles in the world: core, semi-periphery and periphery 
Accordingly, in the visualized network, the countries labeled by them are distributed along three 
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circles from inside to outside, respectively. 

A two-layer hierarchical organization has been constructed by our system, as illustrated in 
Fig. [3]± A reduced blocking model is shown in the first layer by neglecting trivial couplings, in 
which six blocks are obtained. By referring the matrix of the network as presented in Fig. [3J3, 
one can observe that the nodes within the same blocks demonstrate homogeneous row as well as 
column connection profiles. By referring to their couplings, ten isomorphism subgraphs of the 
patterns as defined in Fig. [2] are recognized, which, respectively, are one authority, four communi- 
ties, three hubs and two outliers, as shown in Fig. [3^. Quite a few interesting things can be read 
from these uncovered multiplex structures. Globally, the countries near center tend to have larger 
out-degrees, while those far from center have smaller even zero out-degrees. Specifically, (1) ac- 
cording to the coupling strength from strong to weak, three detected hubs can be ranked into the 
sequence of blocks 4, 1 and 3. The first two hubs consist of the "core" of the trade system except 
for Spain and Denmark; (2) the countries from blocks 3, 2 and 6 consist of the backbone of the 
"semi-periphery"; (3) more than a half of "periphery" countries (10 of 17) have zero out-degrees; 
(4) interestingly, the detected blocks are also geography-related, as illustrated in Fig. ^p. Most 
countries of hubs 4 and 1 are from western Europe expect America, China and Japan; most of hub 
3 are from southeastern Asia; most of the community block 6 are from north or south America; 
most of the outlier block 2 are from eastern Europe; most of the outlier block 5 are from Africa or 
some areas of Asia. 

In the second layer, a macroscopic hub-outlier pattern with strong couplings is recognized. 
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Hub blocks 4 and 1 in the first layer collectively form a more global hub of the whole network 
as the core of the entire trade system; other blocks form a more global outlier of the network 
corresponding to the semi-periphery and periphery of the trade system. This nested hub-outlier 
patterns perhaps give us an evidence about how cash flowed through the world in different levels 
in 1994. Note that arc i —> j denotes that country i imported commodities from country j, which 
also indicates that the spent cash has flowed from i to j. In this way, one can image cash flows 
along these arcs from one country to another. According to the global pattern in the second layer, 
the dominant cash flux will flow from the core countries to themselves with a probability 0.6, and 
to the rest with a probability 0.57. Locally, the blocking model in the first layer shows the backbone 
of the cash flow through the entire world with their respective strength in terms of block couplings, 
as illustrated in Fig. |3ji. 

Mining granular association rules from networks When a network encodes the co-occurrence 
of events, its underlying node- or block-coupling matrices would imply the probabilistic asso- 
ciations among these events in different granular levels, respectively. Formally, we have: node 
association rule (NAR): i j (pij), and block association rule (BAR): B p — )► B q {(j) pq ). A NAR 
means that event j would happen with a probability p^ given event i happens. A BAR means that 
any event of block q would happen with a probability c/) pq given any event of block p happens. Such 
association rules with different granularities can be used in making prediction in a wide range of 
applications, such as online recommender systems. We will demonstrate this idea with a political 
book co-purchasing network compiled by V. Krebs, as given in Fig. |4|a), where nodes represent 
books about US politics sold by the online bookseller Amazon.com, and edges connects pairs of 
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books that are frequently co-purchased, as indicated by the "customers who bought this book also 
bought these other books" feature on Amazon. Moreover, these books have been labeled as "lib- 
eral", "neutral" or "conservative" according to their stated or apparent political alignments based 
on the descriptions and reviews of the books posted on Amazon.com ^2. 

A two-layer hierarchical organization has been detected by our system as shown in Fig. |4j 
The blocking model of the first layer is shown in Fig. |4|b). By matching isomorphism subgraphs in 
its reduced blocking model, nine patterns are recognized, which respectively are five communities 
(blocks 1,2,4,6,7), two cores (blocks 2 and 7), two outliers (blocks 3 and 5) and a bow tie (blocks 
1,2,3). Note that, in indirected networks, the core of a bow tie (block 2 in this case) can be seen 
as the overlapping part of its two wings (blocks 1 and 3 in this case) by neglecting the direction 
of links. In the second layer, a macroscopic community structure is recognized, as shown in Fig. 
|4|c). Interestingly, the left and right communities can be globally labeled as "left-wing" and "right- 
wing" according to the types of the books they contain respectively. Such a global pattern can be 
seen as one macroscopic BAR: the books with common labels would be co-purchased with a great 
chance (about 15%); while, those with different labels are rarely co-purchased (only with a chance 
of 1%). 

When zooming in to both global communities in the second layer, one will obtain 7x7 
mesoscopic BARs encoded by the block-coupling matrix $, as illustrated by the weighted arrow 
lines Fig. |4|b). As an example, we list the BRAs related to the block 2 in a decreasing sequence 
of association strength. B 2 £ 2 (0.60); B 2 £ 3 (0.44); B 2 Si (0.21); B 2 -> £ 4 (0.06); 
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B 2 —> 5 6 (0.04); B 2 —> 5 5 (0); B 2 —> B 7 (0). Such mesoscopic association rules would help 
booksellers adaptively adjust their selling strategies to determine what kinds of stocks they should 
increase or decrease based on the statistics of past sales. For example, if they find the books labeled 
as block 2 are sold well, they may correspondingly increase the order of books labeled as blocks 
1 and 3 besides block 2, while they may simultaneously decrease the order of books labeled as 
blocks 5 or 7. 

More specifically, with the aid of NARs, booksellers would be able to estimate the chance that 
customers will spend their money on a book j if they have already bought book i by referring to i —> 
j(Pij). Such microscopic rules would provide booksellers the suggestions on what specific books 
should be listed according to what sequence in the recommending area of the web page advertising 
a book. For example, for the book "The Price of Loyalty" , the most worth recommended books 
are listed as follows according to the coupling strength to it: Big Lies(0.91); Bushwhacked (0.73); 
Plan of Attack(0. 73); The Politics of Truth (0.7 '3); The Lies of George W.Bush(0.73); American 
Dynasty (0.64); Bushwomen (0.64); The Great Unraveling (0.64); Worse Than Watergate (0.64). 

5 Conclusions 

In this work, we have shown through examples that the structural patterns coexisting in the same 
real-world complex network can be miscellaneous, overlapped or nested, which collaboratively 
shape a heterogeneous hierarchical organization. We have proposed a framework based on the 
concept of granular couplings and the proposed granular blocking model to uncover such multi- 
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plex structures from networks. From the output of patterns, hierarchies and granular couplings 
generated by our approach, one can analyze or even predict some dynamics of networks, which 
are helpful for both theoretical studies and practical applications. 

Moreover, based on the rationale behind this work, we suggest that the evolution of a real- 
world network may be driven by the co-evolution of its structural patterns and its underlying cou- 
plings. On one hand, statistically significant patterns would be gradually highlighted and emer- 
gency shaped by the aggregations of homogeneous individuals in terms of their couplings. On 
the other hand, individuals would adaptively adjust their respective couplings according to the 
currently evolved structural patterns. 
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Support Information 



A Proofs and algorithms 

Proposition 1. For an indirected network, its feedforward-coupling matrix P is equal to its feedback- 
coupling matrix Q. 

Proof 

Pij = P(i-> j\y = k) 

= Ef=i m ik P(i^j\y = k) 

= Ef=i(Ef=i hzik)P(i -+j\y = k) 

= Ef=i EfeLi hzikOkj 

where i —> j denote the event that node V{ couples with Vj, and y = k denote the event that V{ is 
labeled by cluster k\ m ik = 1 if Vi is labeled by cluster k, otherwise m ik = 0. So we have 

p = B 9 ze. 

q i:j = P(i <•— = fc) 
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= Ef=i(Ef=i hzik)P(i +— j\y = k) 



- Ya=1 Ylk=l b il Z lk$i 



kj 



where i <- — j denote the event that node V{ except to be coupled by vj. So we have 

Q = B g ZA. 



If A is symmetry, from the Eq.4 in the article, we have 

2^1 = 1 Z^fe.^O ilk 2^1=1 2^b il7 ^o ilk 

So we have P = Q. □ 



Proposition 2. 



K 



L(N\X,B g ) = J2Y1 ln E n M(^%)/fe a i>fc ( 1Q ) 

1=1 bu^O k=l 

where /(x, y) = — x) l ~ y . 

Proof. Let i> = i denote the event that a node with linkage structure < an, • • • , a in , an, • • • , a ni > 
will be observed in network N. Let y = k denote the event that the cluster label assigned to a node 
is equal to k. Let i —> a j denote the event that node vi link to node Vj or not depending on a^. 
Let i ^aji j denote the event that node Vi will be linked by node vj or not depending on a^. We 
have: 

L(N\X, B g ) = In U? =1 P(v = i) = YT i=l In P(v = i) 
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= E^iInP((t; = t)n(U^ 1 y = fe)) 

= ELi ln Ef=i p ( v = i,y = k) 

= El 1 ^Ek= 1 (P(v = i \y = k)P(y = k)) 

= Er=i ln Ef=i(P(< a n, ■ ■ ■ ,ain,a u , ■■■ ,a ni > \y = k)P(y = k)) 

= Er=i inEf=i(n^ =1 P(i j'b = fe)P(< <-«„ j|y = fc)P(» = fc)) 
= Er=i inEf=i(n^ =1 (^7(i - e kj y-^)(5%(i - s^-*)^) 

= Er=i ln Ef=i( n j=i/(^. Oij)f(Skj, aji)oj k ) 

= Ef=i E 6iI #o ln Ef=i( n "=i/(^i, <kj)f(S kj , aji)u k ) □ 
Proposition 3. 

L K n 

L(JV,Z|X,5 ff ) = Yl Z^E^/^ay) + ln /fe%*)) + ln ^) ( n ) 

/=1 6^0 fe=l j=l 

Proo/ Let denote the cluster label assigned to node i under the given partition Z, we have: 

L(7V,Z|X,£g 

= InnjL 1 P(v = i,y = y(i)) 

= Er=i ln Ef=i ^P(^ = i,y = k) 
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= Eli ln E fc =i m ik P(v = i\y = k)P(y = k) 
= Eli Ef=i HP(v = i\y = k)P(y = k)) mik 

= Eli Ef=i ^ H p ( v = Av = k ) p (y = k )) 

= Eli sL^M^i^a - C'W 1 - 

= Eli Ek=i ^(E^iO 11 /^, ay) + ln/(4 i5 o^)) + lnu; fc ) 
= Eb tl ^ Ek=i z ik(E]=i( l n f(0kj, dij) + lnf(5 kj , aji)) + ha.Uk) + ■ ■ ■ 
+ Eb tL ^ Ef=i ZLk{E"=i( lli f(°kj, dij) + ln/(<*jy, aji )) + lnu k ) 
= Ef=i E 6iI #o Ef=i ^fc(E"=i( ln /(^, Oij) + ln/(4i, aji)) + Inwfc). 

□ 

Notice that, in the proofs of Props [2] and |3} all probabilities such as P(y = k\v = i) and 
P(y = k) are discussed under the conditions of X and B g . To simplify the equations, we omit 
them without losing correctness. 

Proposition 4. 

^ _ 1 ^j=if(hj, aij)f(S kj , aji)u k ^ 

Ei=i hi , , ^ fc=1 I^ =1 f(d kj ,a ij )f(8 kj , aji)u k 
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Proof, let P(y = k\v = i) be the probability that node i belongs to cluster k given X and B g . We 
have: 

lik = P(y = k\b = /, X, B g ) = E^/o ^tb- t P (y = = z ) 
where ^n 1 ^ is the probability of selecting node i from block /. 

According to the Bayesian theorem, we have: 



P(y = k\v = i) 



P(y=k)P(v=i\y=k) 
^ =1 P(y=k)P(v=i\y=k)- 



Based on the proof of Propj2} we have: 

P(y = k)P(v = i\y = k) = U] =1 f(9 kj , a ij )f(5 kj) a ji )uj k . 

So, we have Eq{l2} □ 

As an approximate version of Eq|l2} we have: 

lik = P(y = k\b = l) = P(y = k\v = i, b u ^ 0) 

= P(y = k)P(v = i,b tl ^0\y = k) 

Ef=i P(y = k)P(v = z, b u ^ 0\y = k) 

\i^o^j=if(0kj^ij)f{S kj , aji)uj k 

— k (13) 

Z)fc=i \i^o^j=if(0 kj , aij)f(5 kj , aji)uj k 

where 3&. z ^ denotes that randomly selecting a node from block /. 



That is, instead of averaging all nodes in the block /, the real value of ji k can be approxi- 
mately estimated by a randomly selected node from block I. 
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Correspondingly, an approximate version of the log-likelihood of EqjTUJis given by: 



K 



L(N\X,B g ) = J2 N i(\ l ^Y, U Uf( kj, a ij)f( S k j ,a ji )u; k ) (14) 



1=1 k=l 



where Ni denotes the size of block I. 



The time calculating Eqs|l2]and 10 will be bounded by 0{n 2 K). While, the time calculating 



Eqs[l3] and 14 will be bounded by O(LnK). This will be much efficient for constructing the 



hierarchical organizations of networks. 

Theorem 1. A local optimum ofEq^l^will be guaranteed by recursively calculating Eqs.4 and 5 
in the article. 

Proof. From the Propositioi^2j we have: 

L(N\X,B 9 ) 

= El^P(v = t\X,B 9 ) 

= El 1 ^Ek= 1 P(v = hy = k\x,B 9 ) 

- V n 1tiV K Pin — k\ii — i p \ P(v=i,y=k\X,B g ) 

(by Jensen's inequality) 
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Furthermore, we have: 

G(XW,XW) 

_ v-n v-K „/ _ , | _ • ^ ( S ) o \ i P(v=i,y=k\xW,B g ) 
- 2^i=i 2^k=i r \y rC\V I, A. , £> g ) in P(2/=fc |„ =i;X ( s ) ;B9) 

= Er=i Ef=i P(V = k\v = i, IW,5 S ) In P(„ = i|xW, 5,) 

= Er=i in p(« = i\x(°\ B g ) yLi p (y = k \ v = * (s) ,^)) 

= Er=iln^ = ^ W ,£ fl ) 
= L(iV|X( s ),5 s ). 

Let P(y = k\b = I, X (s \B g ) = 7^, we have: 

G(X,X^) 

= Ef=iE fc ^oEf=i7; fc s) lnP(^ = i,y = k\X,B g ) - Ef=iE^oEf=i7^1nP(y = = 
So, we have: 
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argmaxG(X,XW) 

= argmax(£f=i £ 6 . i#0 E£=i 7/j? h^O = i,y = k\X, B g ) - £f =1 E^o E?=i 7a? lnP (v = 
k\v = i,X( s \B g )) 

= argmax(^f =1 Y.b«^ YlLii^k lnP ( v = hV = k\ x , B g))) 
= arg max E[L(N, \X, B g )] 
= X {s+1 \ 

Recall that, the 6 (s+1) , /\ (s+1) and ft( s+1 ) of X^ s+1 ^ can be computed in terms of ^ by Eq.4 
in the article. So, we have: 

G(X( S+1 \X^) > G(X< s U (s) ) = L(7V|XW,S fl ). 

Recall that L(N\X, B g ) > G(X, X&), we have: 

L(N\X( S+1 \ B g ) > G(X^ S+1 \X^) > G(X( S \X^) = L(N\X^ s \B g ). 

That is to say, the X^ s+1 ^ obtained in the current iteration will be not worse than X^ obtained 
in last iteration. So, we have the theorem. □ 

Proposition 5. In terms of the parameter of X, 0, A, Z and Q, we have: 

$ = QBgZD' 1 , # = ABgZD' 1 (15) 

where D = diag(ntt). 
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Proof. We have 

$pq = Y^ieC q W^pi 

where i e C q denotes node i is in the cluster q with a size N q , and is the probability of selecting 
node i from cluster q. Furthermore, we have: 

Similarly, we have: 

^PQ = i SLl 5pi{BgZ) iq . 

So, we have 

$ = OBgZD- 1 , * = AB g ZD-\ 

□ 

Algorithm 1. Searching the optimal model X given N and B g 
X-GBM(N,B g ) 

01. K = 1; 

02. = LM(N, K, B g ); = - \n(N\X^°\ B g ); 

03. forK=2:n 

04. X« = LM(N, K, B g ); 
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05. = -\n(N\X^,B g ) + 2gK 2 lnK 2 ; 

06. ifLM < L (0) then 

07. X® = X«; L(°) = L«; 
OS. enJj/ 

09. if keeps increasing for predefined steps then 

10. return X^; 

11. endif 

12. endfor 

13. return X (0) ; 

Algorithm 2. Searching the optimal model X given N, K and B 9 
X=LM(N,K,B g ) 

01. randomly initialing = ("fi^LxK with J^k Ilk = 1>' 

02. s <- 1; 

34 



03. repeat until convergence: 

04. compute 6^ and 

05. compute T^; 

06. s <- s + 1; 

07. compute Z from according to Bayesian rule; 

Proposition 6. Let B 9i denotes the blocking model on the i — th layer of the hierarchical organi- 
zation of network N, we have: 

P(X\N, B 91 ,---, BJ ex P(N\X, B 9h )P(Xy» 
Proof. P(X\N,B 9l ,--- ,BJ 

P(X,N,B gi ,-,B 9h ) 
P(N,B 91 ,-,B 9h ) 

<xP(X,N,B gi ,--- ,B gh ) 

= P{N\X,B (n ,... ,B gh )P(X,B gi ,-.. ,B 9h ) 

= P(N\X,B 91 ,-.. ,B 9h )P(X\B 91 ,... ,B 9h ) 

P(B 9h \B 9l ,- ■ ■ , fl^J • • • P(B 92 \B 91 )P(B 91 ) 
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ocP(N\X,B gi ,-.. ,B gh )P(X\B gi ,-.. ,B 9h ) 

Since two nodes from the same block of B a . , will also be in the same block of B a . , we have: 

y%— l yi ' 

P(X|7V, B 9l , • • • , flj oc P(7V|X, S,JP(X|^J = P(7V|X, B 9h )P(X)^ 

□ 

Algorithm 3. Computing the threshold of a blocking model based on $ 
07. sort all entries ofQ into a non-increasing sequence S; 

02. cluster all entries of& by the remarkable gaps of S; 

03. return the biggest entry of the cluster with the minimum 
mean as the threshold; 

As examples, Fig. [5] illustrates, by means of the above algorithm, how to choose reasonable 
thresholds for reducing blocking models of the world trade net and the co-purchased political book 
network as discussed in the article. 

B Additional Examples 

Community detection from a benchmark network We have applied our approach to the football 
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association network, a benchmark widely used for testing the performance of community detection 
algorithms. Fig. [6] gives the experimental results. 

Our approach automatically find out 10 clusters or communities from this network. Fig. [(Ja) 
shows the clustered network, and Fig. [(Jcl) shows the clustered matrix, in which dots denote the 
non-zero entries, and the rows and column corresponding to the same cluster will be put together 
and be separated by the solid lines. Fig. [6}} shows the searching process, in which the cost in 
terms of — L(N\X, Bi) + L^(X) firstly drops down quickly, and then goes up sharply after a 
local minimum corresponding to K = 10. This iterative searching process shapes a well-like 
landscape by penalizing both small likelihood and big coding complexity or smaller priori, and a 
good compromise between them is what we expect. Notice that, this estimated community number 
is a little bit smaller than the real one, in total 12 real associations, among which the teams from 
two small independent associations prefer to play matches with outside teams, (c)-(e) show the 
nested blocking models of this network. A hierarchical organization with three layers have been 
constructed. 

Testing against synthetic networks We have tested our approach against several synthetic net- 
works, and here we give two examples. Fig. [7^ give one synthetic network, in which a 2- 
community pattern and two bipartite patterns coexist together. A two-layer nested organization 
are found out. In the first layer, six clusters with homogeneous row as well as column connection 
distributions are detected; in the second layer, such detected clusters are grouped in pairs accord- 
ing to the node couplings with coarser granularity. This time, a 2-community pattern, a reciprocal 
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bipartite pattern and an unilateral bipartite pattern are emerged. Fig. [7]1 give another synthetic 
network, in which 12 communities are organized in a 3 -layer hierarchical structure according to 
the density of connections. Correspondingly, a three-layer nested blocking model is constructed as 
given in Fig. [7^. 
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Figure 1. An example of multiplex structures consisting of hubs, outliers, and communities 
co-existing in one social network. The network as shown in (a) depicts the co-appearances of 
77 characters in the novel Les Miserables (by Victor Hugo). Nodes denote characters and links 
connect any pair of characters that appear in the same chapter of the novel. This data set is from 
The Stanford GraphBase: A Platform for Combinatorial Computing, edited by D.E.Knuth 17 \ The 
physical connection profiles of nodes are represented in terms of the adjacency matrix as shown 
in (b), in which dots denote "1" entries. In total, six blocks are detected by our method in terms 
of the connection profiles of nodes, as separated by solid lines, so that the nodes within the same 
blocks will demonstrate homogeneous connectivities. To clearly visualize the block organization, 
the matrix has been transformed by putting together the nodes within the same blocks which are 
labeled by "block 1" to "block 6" from the top down. Correspondingly, in the network as shown in 
(a), nodes are also colored according to their block IDs (specifically, the coloring schema is: block 
1-red, block 2-green, block 3-blue, block 4-cyan, block 5-gray and block 6-yellow), and the sizes 
of nodes are proportional to their respective degrees, (c) shows the blocking model of the network, 
in which each block is colored, numbered and placed according to the nodes it contains, and the 
sizes of blocks are proportional to the number of nodes they contain, respectively. The weights 
on the arrow lines globally measure the probabilities that one node from one block will connect to 
another from other blocks, (d) shows the reduced blocking model by cutting the arrow lines with 
trivial weights. In this case, one hub pattern consisting of block 1, one outlier pattern consisting of 
block 5 and four community patterns consisting of blocks 2,3,4 and 6 will readily be recognized 
by referring to the probabilistic linkage among these blocks. 
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Figure 2. The schematic illustrations of multiplex structures, (a)-(g) shows seven structural 
patterns frequently observed in real- world networks, which are represented by networks, matrices 
and blocking models, respectively. In the matrices, shades represent the densities of links. In the 
blocking models, circles denote blocks and solid arrow lines denote block feedforward-couplings. 
From (a)-(g), structural patterns are communities, authorities and outliers, hubs and outliers, a bow 
tie, a multipartite, a bipartite and a hierarchical organization, respectively, (g) shows a two-layer 
hierarchy; two overlapped communities (also can be seen as a hub with two communities) and one 
bipartite are in the first layer; two communities are in the second layer. 

Figure 3. Detecting the nested core-periphery organization and the cash flow patterns from a 
world trade net. The network as shown in (a) encodes the trade relationship among 80 countries. 
Nodes denote countries and arcs denote countries imported commodities from others. The physical 
connection profiles of nodes are represented in terms of its adjacency matrix as shown in (b), in 
which dots denote "1" entries. In total, six blocks are detected as separated by solid lines so that 
the nodes within the same blocks will demonstrate homogeneous row- and column- connection 
distributions. As before, the matrix has been transformed by putting together the nodes within 
the same blocks which labeled by "block 1" to "block 6" from the top down. Correspondingly, 
in the network shown in (a), nodes are colored according to their block IDs (specifically, block 
1 -green, block 2-yellow, block 3-cyan, block 4-red, block 5-gray and block 6-blue), and the sizes 
of nodes are proportional to their respective out-degrees, (c) shows the world map in which 80 
countries are colored by the same coloring schema defined above, (d) shows the detected two- 
layer hierarchical organization, in which each block is colored by the same coloring schema defined 
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above and the sizes of blocks are proportional to the number of nodes they contain, respectively. 
In the first layer, a reduced blocking model is given by cutting the trivial block couplings. The 
arrow lines with weights denote the cash flows from the countries within one block to others. 
The cash flows from two hubs (blocks 4 and 1) are highlighted by thicker arrow lines, of which 
the thickness is proportional to their strength measured by respective block couplings. These 
highlighted cash flows would outline the backbone of the cash circulation in the world trade system. 
A macroscopic hub-outlier pattern on the second layer is detected; together with the microscopic 
hub-outlier patterns on the first layer, the whole trade system would demonstrate a nested core- 
periphery organization, (e) shows the multiplex structures discovered from the reduced blocking 
model on the first layer by a procedure of matching isomorphism subgraphs as defined before. 
These structural patterns are stored in their respective reservoirs labeled as authority, community, 
hub or outlier. 

Figure 4. Mining granular association rules from a co-occurrence net. The network in (a) 
encodes the co-purchased relationship among 105 books, in which the nodes with large, median 
and small sizes are labeled by "liberal", "conservative" and "neutral(including one unlabeled)", 
respectively. A two-layer hierarchical organization is detected. In the first layer, seven blocks 
are detected. As before, the matrix as shown in (b) has been transformed by putting together 
the nodes within the same blocks which labeled by "block 1" to "block 7" from the top down. 
Correspondingly, in the network shown in (a), nodes are colored according to their block IDs 
(specifically, block 1-blue, block 2-red, block 3 -gray, block 4-brown, block 5 -green, block 6-cyan 
and block 7-yellow). In the blocking model as shown in (b), each block is colored by the same 
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coloring schema defined above, the sizes of blocks are proportional to the number of nodes they 
contain, respectively; the arrow lines to be reserved in its reduced blocking model are highlighted 
by thicker arrow lines. The blocking model and its corresponding matrix of the second layer 
are shows in (c), in which a global two-community structure is detected. The granularity of the 
network corresponding to different layers is given by g. 
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Figure 1 : An example of the multiplex structures in a social network. 
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Figure 2: The schematic illustrations of multiplex structures. 
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Figure 3: Detecting the nested core-periphery organization and the cash flow patterns from a world 
trade net. 
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(a) ground layer (g=1) 




(b) first layer (g=15) 




10 20 30 40 50 60 70 80 90100 



(c) second layer (g=52.5) 



Figure 4: Mining granular association rules from a co-occurrence net. 
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Figure 5: The illustrations of calculating thresholds for reducing blocking models. The top: cal- 
culating the threshold for the blocking model of the first layer of the world trade net. The bottom: 
calculating the threshold for the blocking model of the first layer of the co-purchased political 
book network. In both figures, the x denotes the sorted entries and the y denotes the values of 
block couplings in terms of $. These couplings are clustered into three or four groups, as sepa- 
rated by rectangles, by remarkable gaps. The calculated thresholds are the maximum entries of the 
clusters with minimum means pointed out by red solid lines, respectively. 
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Figure 6: Detecting communities from the football association network, (a) The clustered football 
association network; (b) The iterative searching landscape. (cl-c2) The matrix and blocking model 
of the first layer. (dl-d2) The matrix and blocking model of the second layer. (el-e2) The matrix 
and blocking model of the third layer. The coloring schemes used in (a) and (c2), (dl) and (el) are 
same. 
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Figure 7: Testing synthetic networks, (a) The matrix of a synthetic network, in which dots denote 
the non-zero entries, and the six detected clusters are separated by solid lines; (b) The blocking 
model in the first layer, in which blocks are numbered according to the same sequence of the clus- 
ters in the matrix, (c) The blocking model in the second layer, (d) Another synthetic network, 
in which twelve detected clusters are separated by the solid lines, (e) The obtained nested block- 
ing model corresponding to the hierarchical organization with three layers, in which blocks are 
numbered according to the same sequence of the clusters in the matrix. 
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